Установка и подключение необходимых библиотек
#install.packages('gsheet', dep = T)
#install.packages('tidyverse', dep = T)
library(gsheet)
library(tidyverse)
library(lubridate)Квартет Энскомба — четыре набора числовых данных, у которых простые статистические свойства идентичны, но их графики существенно отличаются. Каждый набор состоит из 11 пар чисел. Квартет был составлен в 1973 году английским математиком Ф. Дж. Энскомбом для иллюстрации важности применения графиков для статистического анализа, и влияния выбросов значений на свойства всего набора данных.
Эти данные состоят из четырёх пар x и y с практически равным средним значением (M[xi]=9, M[yi]=7.5) и дисперсией между соответствующими элементами пар (D[xi]=11, D[yi]≈4.13) , а также равным коэффициентом корреляции (cor(xi,yi)=0.816). Модель линейной регрессии, построенная методом МНК для всех вариантов описывается уравнением y=3.00+0.500x [2].
Графики представлены на рисунке ниже , из которого видно, насколько могут различаться данные, описываемые внешне статистически одинаково.
require(stats); require(graphics)
asc <- anscombe
ascx1 <dbl> | x2 <dbl> | x3 <dbl> | x4 <dbl> | y1 <dbl> | y2 <dbl> | y3 <dbl> | y4 <dbl> |
|---|---|---|---|---|---|---|---|
| 10 | 10 | 10 | 8 | 8.04 | 9.14 | 7.46 | 6.58 |
| 8 | 8 | 8 | 8 | 6.95 | 8.14 | 6.77 | 5.76 |
| 13 | 13 | 13 | 8 | 7.58 | 8.74 | 12.74 | 7.71 |
| 9 | 9 | 9 | 8 | 8.81 | 8.77 | 7.11 | 8.84 |
| 11 | 11 | 11 | 8 | 8.33 | 9.26 | 7.81 | 8.47 |
| 14 | 14 | 14 | 8 | 9.96 | 8.10 | 8.84 | 7.04 |
| 6 | 6 | 6 | 8 | 7.24 | 6.13 | 6.08 | 5.25 |
| 4 | 4 | 4 | 19 | 4.26 | 3.10 | 5.39 | 12.50 |
| 12 | 12 | 12 | 8 | 10.84 | 9.13 | 8.15 | 5.56 |
| 7 | 7 | 7 | 8 | 4.82 | 7.26 | 6.42 | 7.91 |
##-- now some "magic" to do the 4 regressions in a loop:
ff <- y ~ x
mods <- setNames(as.list(1:4), paste0("lm", 1:4))
for(i in 1:4) {
ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
## or ff[[2]] <- as.name(paste0("y", i))
## ff[[3]] <- as.name(paste0("x", i))
mods[[i]] <- lmi <- lm(ff, data = anscombe)
print(anova(lmi))
}## Analysis of Variance Table
##
## Response: y1
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 27.510 27.5100 17.99 0.00217 **
## Residuals 9 13.763 1.5292
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Variance Table
##
## Response: y2
## Df Sum Sq Mean Sq F value Pr(>F)
## x2 1 27.500 27.5000 17.966 0.002179 **
## Residuals 9 13.776 1.5307
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Variance Table
##
## Response: y3
## Df Sum Sq Mean Sq F value Pr(>F)
## x3 1 27.470 27.4700 17.972 0.002176 **
## Residuals 9 13.756 1.5285
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Variance Table
##
## Response: y4
## Df Sum Sq Mean Sq F value Pr(>F)
## x4 1 27.490 27.4900 18.003 0.002165 **
## Residuals 9 13.742 1.5269
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## See how close they are (numerically!)
sapply(mods, coef)## lm1 lm2 lm3 lm4
## (Intercept) 3.0000909 3.000909 3.0024545 3.0017273
## x1 0.5000909 0.500000 0.4997273 0.4999091
lapply(mods, function(fm) coef(summary(fm)))## $lm1
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0000909 1.1247468 2.667348 0.025734051
## x1 0.5000909 0.1179055 4.241455 0.002169629
##
## $lm2
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.000909 1.1253024 2.666758 0.025758941
## x2 0.500000 0.1179637 4.238590 0.002178816
##
## $lm3
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0024545 1.1244812 2.670080 0.025619109
## x3 0.4997273 0.1178777 4.239372 0.002176305
##
## $lm4
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0017273 1.1239211 2.670763 0.025590425
## x4 0.4999091 0.1178189 4.243028 0.002164602
## Now, do what you should have done in the first place: PLOTS
op <- par(mfrow = c(2, 2), mar = 0.1+c(4,4,1,1), oma = c(0, 0, 2, 0))
for(i in 1:4) {
ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
plot(ff, data = anscombe, col = "red", pch = 21, bg = "orange", cex = 1.2,
xlim = c(3, 19), ylim = c(3, 13))
abline(mods[[i]], col = "blue")
}
mtext("Anscombe's 4 Regression data sets", outer = TRUE, cex = 1.5)Еще один интересный датасет на эту тему. И опять таки - все статистики очень близки. Но если их визуализировать…
library(datasauRus)
datasaurus_dozendataset <chr> | x <dbl> | y <dbl> | ||
|---|---|---|---|---|
| dino | 55.38460 | 97.17950000 | ||
| dino | 51.53850 | 96.02560000 | ||
| dino | 46.15380 | 94.48720000 | ||
| dino | 42.82050 | 91.41030000 | ||
| dino | 40.76920 | 88.33330000 | ||
| dino | 38.71790 | 84.87180000 | ||
| dino | 35.64100 | 79.87180000 | ||
| dino | 33.07690 | 77.56410000 | ||
| dino | 28.97440 | 74.48720000 | ||
| dino | 26.15380 | 71.41030000 |
if(requireNamespace("dplyr")){
suppressPackageStartupMessages(library(dplyr))
datasaurus_dozen %>%
group_by(dataset) %>%
summarize(
mean_x = mean(x),
mean_y = mean(y),
std_dev_x = sd(x),
std_dev_y = sd(y),
corr_x_y = cor(x, y)
)
}dataset <chr> | mean_x <dbl> | mean_y <dbl> | std_dev_x <dbl> | std_dev_y <dbl> | corr_x_y <dbl> |
|---|---|---|---|---|---|
| away | 54.26610 | 47.83472 | 16.76982 | 26.93974 | -0.06412835 |
| bullseye | 54.26873 | 47.83082 | 16.76924 | 26.93573 | -0.06858639 |
| circle | 54.26732 | 47.83772 | 16.76001 | 26.93004 | -0.06834336 |
| dino | 54.26327 | 47.83225 | 16.76514 | 26.93540 | -0.06447185 |
| dots | 54.26030 | 47.83983 | 16.76774 | 26.93019 | -0.06034144 |
| h_lines | 54.26144 | 47.83025 | 16.76590 | 26.93988 | -0.06171484 |
| high_lines | 54.26881 | 47.83545 | 16.76670 | 26.94000 | -0.06850422 |
| slant_down | 54.26785 | 47.83590 | 16.76676 | 26.93610 | -0.06897974 |
| slant_up | 54.26588 | 47.83150 | 16.76885 | 26.93861 | -0.06860921 |
| star | 54.26734 | 47.83955 | 16.76896 | 26.93027 | -0.06296110 |
if(requireNamespace("ggplot2")){
library(ggplot2)
ggplot(datasaurus_dozen, aes(x=x, y=y, colour=dataset))+
geom_point()+
theme_void()+
theme(legend.position = "none")+
facet_wrap(~dataset, ncol=3)
}Загрузим данные по погоде: Киев, Харьков, Львов, Одесса за период: 01.10.2012 - 27.10.2018
Для тех у кого мало памяти и слабые ноутбуки -
url <- ‘docs.google.com/spreadsheets/d/1rJCuU5NLj7_Z45Hs6vxpjJHvbYDckepKHBTy19K6VyI’
url <- 'docs.google.com/spreadsheets/d/1W6nR5kg44fNaoRi_e_DkcBi-WbX5h9FumtFZX3JkiZM/I'
weather <- read.csv(text=gsheet2text(url, format='csv'), dec =',', stringsAsFactors=FALSE)head(weather,50)dt <int> | date <chr> | city_id <int> | city_name <chr> | lon <dbl> | lat <dbl> | temp <dbl> | temp_min <dbl> | ||
|---|---|---|---|---|---|---|---|---|---|
| 1 | 1349092800 | 01.10.2012 12:00 | 702550 | Lviv | 24.02324 | 49.83826 | 15.00 | 15.00 | |
| 2 | 1349092800 | 01.10.2012 12:00 | 706483 | Kharkiv | 36.25000 | 50.00000 | 18.00 | 13.75 | |
| 3 | 1349096400 | 01.10.2012 13:00 | 703448 | Kiev | 30.51667 | 50.43333 | 17.26 | 17.00 | |
| 4 | 1349096400 | 01.10.2012 13:00 | 702550 | Lviv | 24.02324 | 49.83826 | 15.00 | 15.00 | |
| 5 | 1349096400 | 01.10.2012 13:00 | 698740 | Odessa | 30.73262 | 46.47747 | 22.89 | 22.78 | |
| 6 | 1349182800 | 02.10.2012 13:00 | 706483 | Kharkiv | 36.25000 | 50.00000 | 17.00 | 13.75 | |
| 7 | 1349186400 | 02.10.2012 14:00 | 698740 | Odessa | 30.73262 | 46.47747 | 22.89 | 22.78 | |
| 8 | 1349186400 | 02.10.2012 14:00 | 703448 | Kiev | 30.51667 | 50.43333 | 16.73 | 16.11 | |
| 9 | 1349190000 | 02.10.2012 15:00 | 702550 | Lviv | 24.02324 | 49.83826 | 16.00 | 16.00 | |
| 10 | 1349190000 | 02.10.2012 15:00 | 703448 | Kiev | 30.51667 | 50.43333 | 15.88 | 15.00 |
summary(weather)## dt date city_id
## Min. :1.349e+09 Length:193145 Min. :698740
## 1st Qu.:1.401e+09 Class :character 1st Qu.:698740
## Median :1.455e+09 Mode :character Median :703448
## Mean :1.449e+09 Mean :702802
## 3rd Qu.:1.499e+09 3rd Qu.:706483
## Max. :1.541e+09 Max. :706483
## city_name lon lat temp
## Length:193145 Min. :24.02 Min. :46.48 Min. :-26.00
## Class :character 1st Qu.:30.52 1st Qu.:46.48 1st Qu.: 2.00
## Mode :character Median :30.73 Median :50.00 Median : 10.54
## Mean :30.73 Mean :49.16 Mean : 10.42
## 3rd Qu.:36.25 3rd Qu.:50.43 3rd Qu.: 18.94
## Max. :36.25 Max. :50.43 Max. : 42.22
## temp_min temp_max pressure humidity
## Min. :-26.000 Min. :-26.00 Min. : 962 Min. : 0.00
## 1st Qu.: 2.000 1st Qu.: 2.22 1st Qu.:1010 1st Qu.: 59.00
## Median : 10.000 Median : 11.00 Median :1015 Median : 79.00
## Mean : 9.996 Mean : 10.79 Mean :1015 Mean : 73.15
## 3rd Qu.: 18.000 3rd Qu.: 19.00 3rd Qu.:1021 3rd Qu.: 92.00
## Max. : 42.220 Max. : 47.22 Max. :1054 Max. :100.00
## wind_speed wind_deg clouds_all weather_id
## Min. : 0.000 Min. : 0 Min. : 0.00 Min. :200.0
## 1st Qu.: 2.000 1st Qu.: 80 1st Qu.: 0.00 1st Qu.:800.0
## Median : 3.000 Median :176 Median : 40.00 Median :800.0
## Mean : 3.808 Mean :177 Mean : 41.36 Mean :754.4
## 3rd Qu.: 5.000 3rd Qu.:275 3rd Qu.: 80.00 3rd Qu.:803.0
## Max. :58.000 Max. :360 Max. :100.00 Max. :804.0
## weather_main weather_description
## Length:193145 Length:193145
## Class :character Class :character
## Mode :character Mode :character
##
##
##
Подготовим наши данные
weather$date <- as.Date(weather$date, '%d.%m.%Y %H:%M')
weather$year <- as.factor(year(as.POSIXlt(weather$date, format="%d/%m/%Y")))
oneMonthKiev <- weather %>%
filter(date > '2018-03-01 00:00:00',date < '2018-03-31 00:00:00',city_name == 'Kiev')
oneMonthKievGroup <- weather %>%
filter(date > '2018-03-01 00:00:00',date < '2018-03-31 00:00:00',city_name == 'Kiev') %>%
group_by(date) %>%
summarise(temp = mean(temp))
weatherNum <- select(weather, c(7,10:14))Визуализация одной переменной
plot(weather$temp)plot(oneMonthKiev$date, oneMonthKiev$temp)Параметры xlab и ylab служат для изменения названий осей X и Y соответственно
plot(oneMonthKiev$date, oneMonthKiev$temp,
xlab = 'Дата', ylab = 'Температура', main = 'Температура в Киеве')Параметр type позволяет изменять внешний вид точек на графике. Он принимает одно из следующих текстовых значений:
“p” - точки (points; используется по умолчанию) “l” - линии (lines) “b” - изображаются и точки, и линии (both points and lines) “o” - точки изображаются поверх линий (points over lines) “h” - гистограмма (histogram) “s” - ступенчатая кривая (steps) “n” - данные не отображаются (no points)
plot(oneMonthKievGroup$date, oneMonthKievGroup$temp,
xlab = 'Дата', ylab = 'Температура', main = 'Температура в Киеве',
type = 'l')Эти два параметра контролирут размах значений на каждой из осей графика. По умолчанию они оба принимают значение NULL - в этом случае размах выбирается программой автоматически.
plot(oneMonthKievGroup$date, oneMonthKievGroup$temp,
xlab = 'Дата', ylab = 'Температура', main = 'Температура в Киеве',
type = 'l', ylim = c(1,10))Гистограммы распределения строятся с помощью функции hist(). Чтобы изменить ширину кармана (столбца) гистограммы, необходимо задать параметр breaks =, а цвет задается в параметре col
hist(weather$temp,
xlab = 'Температура', ylab = 'Частота', main = 'Распределение температуры в Украине')hist(weather$temp, breaks = 60,
xlab = 'Температура', ylab = 'Частота', main = 'Распределение температуры в Украине')hist(weather$temp, breaks = 60, col = "lightblue", border = "pink",
xlab = 'Температура', ylab = 'Частота', main = 'Распределение температуры в Украине')hist(weather$temp, breaks = 60, col = "lightblue", border = "pink",
freq = FALSE,
xlab = 'Температура', ylab = 'Частота', main = 'Распределение температуры в Украине')hist(weather$temp, breaks = 30, col = "lightblue", border = "pink",
freq = FALSE,
xlab = 'Температура', ylab = 'Частота', main = 'Распределение температуры в Украине')
lines(density(weather$temp), col = "red", lwd = 2)boxplot(temp ~ city_name, data = weather)boxplot(temp ~ city_name, data = weather,
xlab = "Города",
ylab = "Температура",
main = "Температура по городам Украины",
col = "coral")boxplot(temp ~ city_name, data = weather,
xlab = "Температура",
ylab = "Города",
main = "Температура по городам Украины",
col = c('chartreuse3', 'red' , 'cyan', 'darkorange1'),
horizontal = TRUE)
grid()#install.packages('ggplot2', dep = T)
library(ggplot2)Посмотрим взаимосвязь температуры и давления
ggplot(data = weather) +
geom_point(mapping = aes(x = temp, y = pressure))Тоже самое, но на более разряженных данных за один месяц
ggplot(data = oneMonthKiev) +
geom_point(mapping = aes(x = temp, y = pressure))ggplot(data = weather) +
geom_point(mapping = aes(x = temp, y = pressure,
color = clouds_all))ggplot(data = oneMonthKiev) +
geom_point(mapping = aes(x = temp, y = pressure,
color = weather_main,
#shape = weather_main,
size = wind_speed)) +
labs(
title = "Погода в Киеве",
subtitle = "Март 2018",
caption = "По данным: https://openweathermap.org/",
x = "Температура",
y = "Атмосферное давление",
color = "Погода",
size = "Сила ветра"
)ggplot(data = weather) +
geom_point(mapping = aes(x = temp, y = pressure,
color = weather_main,
size = wind_speed))+
facet_wrap(~weather_main, nrow = 3)ggplot(data = weather) +
geom_point(mapping = aes(x = temp, y = pressure,
color = weather_main))+
facet_grid(city_name~weather_main)ggplot(data = weather) +
geom_point(mapping = aes(x = temp, y = pressure,
color = weather_main))+
facet_grid(.~city_name)ggplot(data = weather,aes(x = temp, y = pressure)) +
geom_smooth()ggplot(data = weather) +
geom_smooth(mapping = aes(x = temp, y = pressure,
linetype = city_name,
color = city_name))ggplot(data = oneMonthKiev) +
geom_point(mapping = aes(x = temp, y = pressure)) +
geom_smooth(mapping = aes(x = temp, y = pressure))ggplot(data = oneMonthKiev, mapping = aes(x = temp, y = pressure)) +
geom_point() +
geom_smooth()ggplot(data = oneMonthKiev, mapping = aes(x = temp, y = pressure)) +
geom_point(mapping = aes(color = weather_main)) +
geom_smooth()oneMonthKiev %>%
ggplot(aes(x = temp, y = pressure)) +
geom_point(mapping = aes(color = weather_main)) +
geom_smooth(data = filter(oneMonthKiev, weather_main == 'Rain'),
se = F
)ggplot(data = weather, mapping = aes(weather_main)) +
geom_bar()ggplot(data = weather, mapping = aes(weather_main)) +
stat_count()ggplot(data = weather) +
geom_bar(mapping = aes(x = weather_main, y = ..prop.., group = 1))ggplot(data = weather, mapping = aes(temp)) +
geom_histogram(binwidth = 1) +
theme_bw()ggplot(data = weather, mapping = aes(weather$temp)) +
stat_bin(binwidth = 1)ggplot(data = weather) +
geom_bar(mapping = aes(weather_main, color = weather_main))ggplot(data = weather) +
geom_bar(mapping = aes(weather_main, fill = weather_main))ggplot(data = weather) +
geom_bar(mapping = aes(weather_main, fill = city_name))ggplot(data = weather) +
geom_bar(mapping = aes(weather_main, fill = city_name),
position = 'fill'
)ggplot(data = weather) +
geom_bar(mapping = aes(weather_main, fill = city_name),
position = 'dodge'
)ggplot(data = oneMonthKiev) +
geom_point(mapping = aes(x = temp, y = pressure)
)ggplot(data = oneMonthKiev) +
geom_point(mapping = aes(x = temp, y = pressure),
position = 'jitter'
)ggplot(data = weather) +
geom_bar(mapping = aes(weather_main, fill = city_name)) +
coord_flip()ggplot(data = weather,mapping = aes(x = city_name,y = temp, fill = city_name)) +
geom_boxplot() +
coord_flip()ggplot(data = weather) +
geom_bar(mapping = aes(weather_main, fill = city_name)) +
coord_polar()Очень удобная библиотека для изучения корреляции между переменными.
library(corrplot)
M <-cor(weatherNum)
corrplot(M, method = "circle", tl.cex=0.6, tl.srt = 45, tl.col = "black", type= "upper", order="hclust")#install.packages('plotly', dependencies = T)
library(plotly)plot_ly(weather, y = ~temp,color = ~year, type = "box")weather %>%
plot_ly(
x = ~year,
y = ~temp,
split = ~year,
type = 'violin',
box = list(
visible = T
),
meanline = list(
visible = T
)
) %>%
layout(
xaxis = list(
title = "Year"
),
yaxis = list(
title = "Temp",
zeroline = F
)
)plot_ly(oneMonthKiev, x = ~temp, y = ~pressure, z = ~humidity, color = ~weather_main) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = 'Temp'),
yaxis = list(title = 'Pressure'),
zaxis = list(title = 'Humidity')))gp <- ggplot(data = weather,mapping = aes(x = temp, y = pressure,
linetype = city_name,
color = city_name)) +
geom_smooth()
ggplotly(gp)Это и не только, на нашем, совместном с Data Science UA, курсе по визуализации данных
Data Science Visualization Course: 2.0.
Приглашаем Вас!